%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Monte-carlo wavefunction simulation for a three-level atom [J. Dalibard et al., PRL68, 580(1992)].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function MonteCarlo

Constant;



G12=23e6;																					%decays	P1/2 - S1/2

G14=.5e6;																					%			P1/2 - D5/2

G24=1e-6;																					%			D5/2 - S1/2



o12=	5*G12;																				%Rabi frequencies

o14= 40*G14;

d12=-.5*G12;																				%detunings

d14=	    0;

dt =   1e-9;

t  =[0:1000]*dt;





% Bloch equation



rho=[1 0 0 0 0 0 0 0 0]';

for jj=1:length(t)

	o12_=o12 *((200<=jj)&(jj<400) || (600<=jj)&(jj<800));

	o14_=o14 *((jj <200) || (400<=jj)&(jj<600) || (800<=jj));

	rho=MonteCarlo_Bloch(dt,rho,G12,G14,G24,o12_,o14_,d12,d14);

	rh1(jj)=rho(1);				rh2(jj)=rho(2);			rh3(jj)=rho(3);

end





% Schroedinger equation



nn1=zeros(size(t));				nn2=nn1;						nn3=nn1;

for j=1:100

   Y =[1 ; 0 ; 0];

   for jj=1:length(t)

	o12_=o12 *((200<=jj)&(jj<400) || (600<=jj)&(jj<800));

	o14_=o14 *((jj <200) || (400<=jj)&(jj<600) || (800<=jj));

		Y=MonteCarlo_Mcwf(dt,Y,G12,G14,G24,o12_,o14_,d12,d14);

		n1(jj)=Y(1)'*Y(1);		n2(jj)=Y(2)'*Y(2);		n3(jj)=Y(3)'*Y(3);

	end

   nn1=nn1+n1;						nn2=nn2+n2;					nn3=nn3+n3;

	MonteCarlo_Plot(j,t,rh1,rh2,rh3,n1,n2,n3,nn1,nn2,nn3);

end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Subroutine Bloch equation



function rho = MonteCarlo_Bloch(dt,rho,G12,G14,G24,o12,o14,d12,d14)



I  =i/2;

g12=G12+G24    ;																			%decays of the coherences

g14=        G14;

g24=G12+G24+G14;

d24=d12-d14;

L12=g12/2+i*d12;

L14=g14/2+i*d14;

L24=g24/2+i*d24;



%11			22			44			12			21			14			41			24			42

M=[0			G12		G14		I*o12		-I*o12	I*o14		-I*o14	0			0

	0			-G12-G24 0			-I*o12	I*o12		0			0			0			0

	0			G24		-G14		0			0			-I*o14	I*o14		0			0 

	I*o12		-I*o12	0			-L12		0			0			0			0			-I*o14

	-I*o12	I*o12		0			0			-L12'   	0			0			I*o14		0

	I*o14		0			-I*o14	0			0			-L14		0			-I*o12	0		

   -I*o14	0			I*o14		0			0			0			-L14'		0			I*o12

	0			0			0			0			I*o14		-I*o12	0			-L24		0

	0			0			0			-I*o14	0			0  	 	I*o12		0			-L24'		];



rho=expm(M*dt)*rho;





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Subroutine Schroedinger equation



function Y = MonteCarlo_Mcwf(dt,Y,G12,G14,G24,o12,o14,d12,d14)



H=[0   		o12/2 			o14/2

	o12/2		-d12-i*G12/1	0

	o14/2		0					-d14-i*G14/1];



Y=expm(-i*H*dt)*Y;		%Y=(eye(2)-i*H*dt)*Y;

Z=rand;

if Y'*Y>Z

	Y=Y/(Y'*Y);

else

	Y=[1 ; 0 ; 0];

end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Subroutine plot



function Y = MonteCarlo_Plot(j,t,rh1,rh2,rh3,n1,n2,n3,nn1,nn2,nn3)



figure(1);												clf;

subplot(311),plot(1e9*t,rh3,'k',1e9*t,n3,'r',1e9*t,nn3/j,'g')

fill([200 400 400 200],[0 0 1 1],'y');			hold on;

fill([600 800 800 600],[0 0 1 1],'y');

plot([200 200],[0 1],'y');							plot([400 400],[0 1],'y');

plot([600 600],[0 1],'y');							plot([800 800],[0 1],'y');

plot(1e9*t,rh3,'k',1e9*t,n3,'r',1e9*t,nn3/j,'g')

hold off;

set(gca,'xLim',[0 max(t)*1e9],'yLim',[0 1],'xTickLabel',' ','FontSize',16,'LineWidth',1.2)

ylabel('\fontsize{16}{\itD}_{5/2}');

text( 50,1.1,'\fontsize{16}{\itS - D}');

text(250,1.1,'\fontsize{16}{\itS - P}');

text(450,1.1,'\fontsize{16}{\itS - D}');

text(650,1.1,'\fontsize{16}{\itS - P}');

text(850,1.1,'\fontsize{16}{\itS - D}');



subplot(312),

fill([200 400 400 200],[0 0 1 1],'y');			hold on;

fill([600 800 800 600],[0 0 1 1],'y');

plot([200 200],[0 1],'y');							plot([400 400],[0 1],'y');

plot([600 600],[0 1],'y');							plot([800 800],[0 1],'y');

plot(1e9*t,rh1,'k',1e9*t,n1,'r',1e9*t,nn1/j,'g')

hold off;

set(gca,'xLim',[0 max(t)*1e9],'yLim',[0 1],'xTickLabel',' ','FontSize',16,'LineWidth',1.2)

ylabel('\fontsize{16}{\itS}_{1/2}')



subplot(313),plot(1e9*t,rh2,'k',1e9*t,n2,'r',1e9*t,nn2/j,'g')

fill([200 400 400 200],[0 0 1 1],'y');			hold on;

fill([600 800 800 600],[0 0 1 1],'y');

plot([200 200],[0 1],'y');							plot([400 400],[0 1],'y');

plot([600 600],[0 1],'y');							plot([800 800],[0 1],'y');

plot(1e9*t,rh2,'k',1e9*t,n2,'r',1e9*t,nn2/j,'g')

hold off;

set(gca,'xLim',[0 max(t)*1e9],'yLim',[0 1],'FontSize',16,'LineWidth',1.2)

xlabel('\fontsize{16}{\itt}   (ns)');

ylabel('\fontsize{16}{\itP}_{1/2}')

j

if j>80

	pause;

end



